---
title: "results of list experiment and balance check"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)

library(tidyverse)
library(haven)
library(stargazer)
library(BalanceR)
library(ggpubr)
library(cobalt)

dat <- read_dta("data/Iraq2023_ver1.dta")

dat <- dat %>% 
  mutate(Gender = sex,
         Age = age,
         Education = education,
         Income = income,
         na.rm = TRUE)
dat <- dat %>% 
  mutate(Distrust_gov = (q7_1_1 + q7_1_2 + q7_1_3 + q7_1_7) / 4,
         trust_gov = (trust_president + trust_PM + trust_parliament + trust_party) / 4) %>% 
  mutate(Distrust_gov2 = case_when(Distrust_gov <= 3 ~ "low",
                                   Distrust_gov >= 3.25 ~ "high"),
         trust_gov2 = case_when(trust_gov <= 3 ~ 0,
                                trust_gov >= 3.25 ~ 1)) %>%
  mutate(Minorities = others)

```

# data summary
```{r}
stargazer(as.data.frame(dat), type = "text")
```

# simple results
```{r}
prop.table(table(dat$violent_control))
prop.table(table(dat$violent_treatment))
prop.table(table(dat$violent_experiment))
```

# balance check
```{r}
balance_violence <- BalanceR(data = dat, group = list_violent,
                         cov = c(Gender, Age, Education, Income, 
                                 Minorities, support_PMU_parties, trust_gov2,
                                 reconciliation, anti_corruption))
print(balance_violence, digits = 3)
summary(balance_violence, digits = 5)
plot(balance_violence, point.size = 5, text.size = 18)

balance_violence2 <- BalanceR(data = dat, group = list_violent,
                         cov = c(Gender, Age, Education, Income))
print(balance_violence2, digits = 3)
summary(balance_violence2, digits = 5)
plot(balance_violence2, point.size = 5, text.size = 18)
```

# t-test
```{r}
# result
mean(dat$violent_treatment, na.rm = TRUE) - mean(dat$violent_control, na.rm = TRUE)
mean(dat$PMU_treatment, na.rm = TRUE) - mean(dat$PMU_control, na.rm = TRUE)

t.test(dat$violent_control, dat$violent_treatment)
t.test(dat$PMU_control, dat$PMU_treatment)

# t-test gender
t.test(Gender ~ list_violent, data = dat)
t.test(Minorities ~ list_violent, data = dat)

# SMO cobalt package
colSums(is.na(dat[, c("list_violent", "Age", "Gender")]))
dat_clean <- na.omit(dat[, c("list_violent", "Age", "Gender")])

bal.tab(list_violent ~ Gender + Age, data = dat_clean)
love.plot(bal.tab(list_violent ~ Gender + Age + Education + Income +
                    Minorities + support_PMU_parties + trust_gov2 +
                    reconciliation + anti_corruption, 
                  data = dat), stat = "mean.diffs")

# regression if treatment is not significant, balanceddat_clean# regression if treatment is not significant, balanced
summary(lm(Gender ~ list_violent, data = dat))
```

# list results
```{r}
list_pmu <- dat %>% 
  summarise(control_mean = mean(PMU_control, na.rm = TRUE), control_sd = sd(PMU_control, na.rm = TRUE),
            treatment_mean = mean(PMU_treatment, na.rm = TRUE), treatment_sd = sd(PMU_treatment, na.rm = TRUE))

sd(dat$PMU_control, na.rm = TRUE) / sqrt(length(dat$PMU_control))
sd(dat$PMU_treatment, na.rm = TRUE) / sqrt(length(dat$PMU_treatment))

p1 <- data.frame(
  Experiment = c("control", "treatment"),
  Mean = c(1.93544, 2.400259),
  SE = c(0.02004044, 0.02671396)
)
p2 <- ggplot(p1, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of PMU") +
  theme_bw()


list_vio <- dat %>% 
  summarise(control_mean = mean(violent_control, na.rm = TRUE), control_sd = sd(violent_control, na.rm = TRUE),
            treatment_mean = mean(violent_treatment, na.rm = TRUE), treatment_sd = sd(violent_treatment, na.rm = TRUE))

sd(dat$violent_control, na.rm = TRUE) / sqrt(length(dat$violent_control))
sd(dat$violent_treatment, na.rm = TRUE) / sqrt(length(dat$violent_treatment))


p3 <- data.frame(
  Experiment = c("control", "treatment"),
  Mean = c(1.873826, 2.2),
  SE = c(0.01766219, 0.02520999)
)
p4 <- ggplot(p3, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of election violent") +
  theme_bw()

p4
p2

ggarrange(p4, p2,
          nrow = 1, ncol = 2, align = "hv")


# new plot
list <- dat |> 
  summarise(control_mean = mean(violent_control, na.rm = TRUE), control_SE = sd(violent_control, na.rm = TRUE) / sqrt(length(violent_control)),
            treatment_mean = mean(violent_treatment, na.rm = TRUE), treatment_SE = sd(violent_treatment, na.rm = TRUE) / sqrt(length(violent_treatment)))

df1 <- data.frame(
  Experiment = c("Control", "Treatment"),
  Mean = c(list$control_mean, list$treatment_mean),
  SE = c(list$control_SE, list$treatment_SE)
)

p5 <- ggplot(df1, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of electoral violence") +
  theme_bw()
p5
ggsave(filename = "result/list_result.png",
       plot = p5,
       width = 4,
       height = 4,
       units = "in", 
       dpi = 600)

```


# policy
```{r}
barplot(table(dat$shia, dat$policy),
        main = "Most important policy for stability",
        names.arg = c("Nationalism", "Patriotism", "Civilization", "Reconciliation", "Democracy", "Justice"),
        beside = TRUE,
        ylim = c(0, 300))
legend("topright", col = c("black", "gray"), pch = 15, legend = c("Shia","others"))   
```

